Libraries:
# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')
Source sapling data:
source("Scripts/SA_Import.R")
Pivot wider to create dataframe where each row is for one plot and has total details for each species group
sa_merge2 <- sapling_merge %>%
pivot_wider(names_from = Species_Groups, values_from = Total_Tally)
Import time since data and add it to the PIRI dataset
time_since <- read_csv("CleanData/Treat_Year_Data.csv")
sa_merge3 <- merge(sa_merge2, time_since, by = 'Site')
#log transform time from treatment data
sa_merge3$l.TFT <- log(sa_merge3$Time_from_Treat)
Run the ‘Add_BA’ script and merge with dataset:
source("Scripts/Add_BA.R")
# merge with sapling dataset -------------------
sa_merge4 <- merge(sa_merge3, prism_BA, by = 'Plot_No')
Run ‘Ground_Data.R’ script and add it to sapling dataset:
source("Scripts/Ground_Data.R")
# merge with sapling dataset -------------------
sa_merge5 <- merge(sa_merge4, ground3, by = 'Plot_No')
Source and add veg data:
source("Scripts/Veg_Data.R")
# merge with sapling dataset -------------------
sa.all <- merge(sa_merge5, veg3, by = "Plot_No")
rm(sa_merge3,
sa_merge4,
sa_merge5)
Sapling count data is taken in 25m2 plots; basal area is measured in hectares; veg and soil data is taken a 1m2 plots.
Sapling data will be convered into 1m2 plots in order to compare across and reduce the amount of scales of data collection, reducing it to two: 1m2 plots and per hectare observations
sa.all$PIRI.1m <- sa.all$PIRI/25
sa.all$SO.1m <- sa.all$Shrub_Oak/25
sa.all$Other.1m <- sa.all$Other/25
Create log transformed categories for newly added variables, then select for just the desired variables:
sa.all$l.PIRI1 <- log(sa.all$PIRI.1m + 1)
sa.all$l.SO1 <- log(sa.all$SO.1m + 1)
sa.all$l.other1 <- log(sa.all$Other.1m + 1)
sa.all2 <- sa.all %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI1, Shrub_Oak, l.SO1, Other, l.other1, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>%
arrange(Treat_Type)
sa.all3 <- sa.all
sa.all3$l.PIRI <- log(sa.all3$PIRI + 1)
sa.all3$l.SO <- log(sa.all3$Shrub_Oak + 1)
sa.all3$l.other <- log(sa.all3$Other + 1)
sa.all3 <- sa.all3 %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI, Shrub_Oak, l.SO, Other, l.other, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>%
arrange(Treat_Type)
Select just for numerical vs log and then look at paired plots:
#not log transformed (1m2)
sa.num <- sa.all2 %>%
select(PIRI.1m, S0.1m, Other.1m, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)
ggpairs(sa.num)
ggpairs(sa.num, aes(color = Treat_Type))
#log transformed (1m2)
sa.numl <- sa.all2 %>%
select(l.PIRI1, l.SO1, l.other1, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)
ggpairs(sa.numl)
ggpairs(sa.numl, aes(color = Treat_Type))
rm(sa.num,
sa.numl)
Can see the correlation coefficients for linear (Pearsons) relationships. None of them appear very strong, except for ones that are analogs (avg LD vs mineral soil exposure; ba/ha vs piri ba/ha)
No significant correlations obvious from plots
#not log transformed (25m2)
sa3.num <- sa.all3 %>%
select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)
ggpairs(sa3.num)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggpairs(sa3.num, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#log transformed (25m2)
sa3.numl <- sa.all3 %>%
select(l.PIRI, l.SO, l.other, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)
ggpairs(sa3.numl)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggpairs(sa3.numl, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
rm(sa3.num,
sa3.numl)
Again no strong correlations
Going to try some scatterplots:
#PIRI & SO
ggplot(sa.all3, aes(x=l.PIRI, y = l.SO))+
geom_point()+
geom_smooth(method = lm)+
facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'
# negative for control, positive for fallrx, flat for other
#PIRI and other
ggplot(sa.all3, aes(x=l.PIRI, y = l.other, fill = Treat_Type))+
geom_point()+
geom_smooth(method = lm) +
facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'
# negative for control & fall rx, weak negative for springrx and mow, positive for harvest
#PIRI and veg
ggplot(sa.all3, aes(x=l.PIRI, y = l.Veg_Total, fill = Treat_Type))+
geom_point()+
geom_smooth(method = lm) +
facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'
# positive for harvest, weak negative for springrx, negative for mowrx, fallrx, and control
#relationships varying by treatment type may be a reason why the model are having a hard time in DHARMa
Going to start with data transformed to 1m plot and no treatment type
sa.m1 <- glmmTMB(PIRI ~ l.SO1 + l.other1 + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all2,
family = poisson)
AIC(sa.m1) #213.7
## [1] 213.6906
sa.m1_sr <- simulateResiduals(sa.m1, n = 1000, plot = T)
testZeroInflation(sa.m1_sr) #zero inflated
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0534, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
sa.m1a <- glmmTMB(PIRI ~ l.SO1 + l.other1 + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all2,
ziformula = ~1,
family = poisson)
AIC(sa.m1a) #215.7
## [1] 215.6909
sa.m1a_sr <- simulateResiduals(sa.m1a, n = 1000, plot = T) #passes
testZeroInflation(sa.m1a_sr) #still zero inflated, lets test on data set without 1m transformations - sa.all3
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0539, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
# Untransformed sapling counts
sa.m1b <- glmmTMB(PIRI ~ l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m1b) #215.9
## [1] 215.8987
sa.m1b_sr <- simulateResiduals(sa.m1b, n = 1000, plot = T) #fails ks
testZeroInflation(sa.m1b_sr) #still very zero inflated
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.056, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
#going to test models with untransformed counts and see if zero inflation issues resolve as variables decrease?
# test piri ba
sa.m2 <- glmmTMB(PIRI ~ l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m2) #217.1
## [1] 217.0877
lrtest(sa.m1b, sa.m2) # p = 0.07
## Likelihood ratio test
##
## Model 1: PIRI ~ l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -96.949
## 2 10 -98.544 -1 3.189 0.07413 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test ba
sa.m3 <- glmmTMB(PIRI ~ l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m3) #215.1
## [1] 215.1096
#test mineral soil
sa.m4 <- glmmTMB(PIRI ~ l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m4) #213.2
## [1] 213.2266
lrtest(sa.m3, sa.m4) # p = 0.7
## Likelihood ratio test
##
## Model 1: PIRI ~ l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: PIRI ~ l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -98.555
## 2 8 -98.613 -1 0.117 0.7323
#test litter depth
sa.m5 <- glmmTMB(PIRI ~ l.SO + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m5) #211.2
## [1] 211.249
#test veg cover
sa.m6 <- glmmTMB(PIRI ~ l.SO + l.other + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m6) #209.4
## [1] 209.3909
lrtest(sa.m5, sa.m6) #p = 0.7
## Likelihood ratio test
##
## Model 1: PIRI ~ l.SO + l.other + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 7 -98.624
## 2 6 -98.695 -1 0.1419 0.7064
#test other
sa.m7 <- glmmTMB(PIRI ~ l.SO + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m7) #253.7
## [1] 253.7722
lrtest(sa.m6, sa.m7) #p less than 0.001
## Likelihood ratio test
##
## Model 1: PIRI ~ l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ l.SO + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 6 -98.695
## 2 5 -121.886 -1 46.381 0.000000000009734 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#test shrub oak
sa.m8 <- glmmTMB(PIRI ~ l.other + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m8) #249.3
## [1] 249.3442
lrtest(sa.m6, sa.m8) #p less than 0.001
## Likelihood ratio test
##
## Model 1: PIRI ~ l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ l.other + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 6 -98.695
## 2 5 -119.672 -1 41.953 0.00000000009348 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#sa.m6 looks best, lets test model
sa.m6_sr <- simulateResiduals(sa.m6, n = 1000, plot = TRUE) #doesn't pass
testZeroInflation(sa.m6_sr) #still fails
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0842, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
#try model with treat type and see if it passes?
sa.m9 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m9) #256.3
## [1] 256.386
lrtest(sa.m6, sa.m9) #p less than 0.001
## Likelihood ratio test
##
## Model 1: PIRI ~ l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 6 -98.695
## 2 10 -118.193 4 38.995 0.00000006983 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sa.m9_sr <- simulateResiduals(sa.m9, n = 1000, plot = TRUE) #passes
testZeroInflation(sa.m9_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.9994, p-value = 0.97
## alternative hypothesis: two.sided
testResiduals(sa.m9_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.031089, p-value = 0.7216
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.2775, p-value = 0.42
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.54
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000662650602409639 )
## 0.002008032
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.031089, p-value = 0.7216
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.2775, p-value = 0.42
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.54
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000662650602409639 )
## 0.002008032
#this model looks great. based on modeling below that starts with TT, I'm going to add veg back into this model and run again
sa.m10 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m10) #218.5
## [1] 218.4779
lrtest(sa.m10, sa.m9) #weird, is now important
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -98.239
## 2 10 -118.193 -1 39.908 0.0000000002662 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sa.m10_sr <- simulateResiduals(sa.m10, n = 1000, plot = TRUE) #now fails
testResiduals(sa.m10_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.083632, p-value = 0.001886
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0000000000023763, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.86
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.005070281
## sample estimates:
## outlier frequency (expected: 0.00160642570281124 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.083632, p-value = 0.001886
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0000000000023763, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.86
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.005070281
## sample estimates:
## outlier frequency (expected: 0.00160642570281124 )
## 0
Model has much better fit and zero inflation issues are resolved when treatment type variable is included. Time to do elimination again but with TT
sa.m1b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m1b) #223.4
## [1] 223.4014
# test piri ba
sa.m2 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m2) #223.9
## [1] 223.9656
lrtest(sa.m1b, sa.m2) # p = 0.1
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 15 -96.701
## 2 14 -97.983 -1 2.5642 0.1093
# test ba
sa.m3 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m3) #222.3
## [1] 222.2865
#test mineral soil
sa.m4 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m4) #220.5
## [1] 220.463
lrtest(sa.m3, sa.m4) # p = 0.7
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 13 -98.143
## 2 12 -98.232 -1 0.1765 0.6744
#test litter depth
sa.m5 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m5) #218.5
## [1] 218.4779
#test veg cover
sa.m6 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m6) #256.4
## [1] 256.386
lrtest(sa.m5, sa.m6) #p less than 0.001
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -98.239
## 2 10 -118.193 -1 39.908 0.0000000002662 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#test other
sa.m7 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m7) #217.6
## [1] 217.5681
lrtest(sa.m5, sa.m7) #p = 0.3
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -98.239
## 2 10 -98.784 -1 1.0902 0.2964
#test shrub oak
sa.m8 <- glmmTMB(PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.m8) #215.6
## [1] 215.6341
lrtest(sa.m7, sa.m8) #p = .8
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -98.784
## 2 9 -98.817 -1 0.0661 0.7971
#so with treatment type included, veg cover is important, but other species of saplings are not
#check model
sa.m8_sr <- simulateResiduals(sa.m8, n = 1000, plot = TRUE) #fails
testZeroInflation(sa.m8_sr) #fails
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.088, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
#try model with so and other saplings?
sa.m5_sr <- simulateResiduals(sa.m5, n = 1000, plot = TRUE)
testZeroInflation(sa.m5_sr) #both still fail ...
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0802, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
I have a very bad headache, so I’m going to stop now. Im so delicate. Anyway, veg seems to messing some stuff up, even if AIC says its important. Take a look at scatter plot/lm of piri and veg for saplings. Then play around with models that have and don’t have it. Could try assinging something other than 1 for the zero inflation…
sa.m10 <- glmmTMB(PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = poisson)
AIC(sa.m10) #221.7
## [1] 221.6967
sa.m10_sr <- simulateResiduals(sa.m10, n = 1000, plot = TRUE)
testZeroInflation(sa.m10_sr) #fails
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0632, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
sa.m11 <- glmmTMB(PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Treat_Type,
family = poisson)
AIC(sa.m11) #223.6
## [1] 223.6341
sa.m11_sr <- simulateResiduals(sa.m11, n = 1000, plot=T)
testZeroInflation(sa.m11_sr) #fails
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0874, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
sa.m12 <- glmmTMB(PIRI ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = poisson)
AIC(sa.m12) #213.9 w/1, 250.5 with region
## [1] 250.5171
sa.m12_sr <- simulateResiduals(sa.m12, n = 1000, plot = TRUE) #fails with 1 as ZI; passes with region as zi
testZeroInflation(sa.m12_sr) #fails with 1 as ZI; passes with region as zi
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99973, p-value = 1
## alternative hypothesis: two.sided
Ok, this is the best fit model so far
sa.m12 <- glmmTMB(PIRI ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = poisson)
AIC(sa.m12) #213.9 w/1, 250.5 with region
## [1] 250.5171
#I'd like to try adding so, other, and veg back into this model, checking AIC and model fit
sa.m13 <- glmmTMB(PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = poisson)
AIC(sa.m13) #221.7
## [1] 221.6967
sa.m13_sr <- simulateResiduals(sa.m13, n = 1000, plot = TRUE) #fails
testZeroInflation(sa.m13_sr) #fails
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0632, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
sa.m14 <- glmmTMB(PIRI ~ Treat_Type + l.SO + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = poisson)
AIC(sa.m14)
## [1] 250.9238
lrtest(sa.m12, sa.m14) #p = 0.2
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 12 -113.26
## 2 13 -112.46 1 1.5933 0.2069
sa.m15 <- glmmTMB(PIRI ~ Treat_Type + l.other + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
#won't converge w/ region
AIC(sa.m15)
## [1] 255.3798
sa.m15_sr <- simulateResiduals(sa.m15, n = 1000, plot = TRUE)
testZeroInflation(sa.m15_sr)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99931, p-value = 0.988
## alternative hypothesis: two.sided
sa.f1 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.f1) #256.4
## [1] 256.386
sa.f2 <- glmmTMB(PIRI ~ Treat_Type + l.other + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
#won't converge w/ region
AIC(sa.f2) #255.4
## [1] 255.3798
sa.f2_sr <- simulateResiduals(sa.f2, n = 1000, plot = TRUE) #passes
testZeroInflation(sa.f2_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99931, p-value = 0.988
## alternative hypothesis: two.sided
testResiduals(sa.f2_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.054361, p-value = 0.1054
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0628, p-value = 0.588
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.52
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000642570281124498 )
## 0.002008032
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.054361, p-value = 0.1054
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0628, p-value = 0.588
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.52
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000642570281124498 )
## 0.002008032
lrtest(sa.f1, sa.f2) # p = 0.3 - can drop SO
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -118.19
## 2 9 -118.69 -1 0.9938 0.3188
sa.f3 <- glmmTMB(PIRI ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~1,
family = poisson)
AIC(sa.f3) #213.9
## [1] 213.9437
sa.f3_sr <- simulateResiduals(sa.f3, n = 1000, plot = TRUE) #fails
testZeroInflation(sa.f3_sr) #fails
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.09, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
lrtest(sa.f2, sa.f3) #p is less than 0.001, but the AIC falls so much, weird
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -118.690
## 2 8 -98.972 -1 39.436 0.000000000339 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sa.f4 <- glmmTMB(PIRI ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = poisson)
AIC(sa.f4) #250.5 with region
## [1] 250.5171
sa.f4_sr <- simulateResiduals(sa.f4, n = 1000, plot = TRUE) #passes
testZeroInflation(sa.f4_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99973, p-value = 1
## alternative hypothesis: two.sided
testResiduals(sa.f4_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.032779, p-value = 0.6584
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.88128, p-value = 0.89
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00100401606425703 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.032779, p-value = 0.6584
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.88128, p-value = 0.89
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00100401606425703 )
## 0
AIC(sa.f4, sa.f2)
## df AIC
## sa.f4 12 250.5171
## sa.f2 9 255.3798
#aic is lower for sa.f4 & more df - therefore better model? maybe graph both to see what is going on?
Plot these babies
library(sjPlot)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
library(sjlabelled)
##
## Attaching package: 'sjlabelled'
## The following object is masked from 'package:forcats':
##
## as_factor
## The following object is masked from 'package:ggplot2':
##
## as_label
## The following object is masked from 'package:dplyr':
##
## as_label
library(sjmisc)
##
## Attaching package: 'sjmisc'
## The following object is masked from 'package:purrr':
##
## is_empty
## The following object is masked from 'package:tidyr':
##
## replace_na
## The following object is masked from 'package:tibble':
##
## add_case
set_theme(base = theme_classic(),
theme.font = 'serif',
axis.title.size = 1.5,
axis.textsize.x = 1.5,
axis.textsize.y = 1.5,
title.size = 2.5,
title.align = "center",
legend.pos = "right",
legend.size = 1.5,
legend.title.size = 1.5,
#legend.bordercol = "black",
legend.item.size = .75)
plot_model(sa.f4,
type = "pred",
terms = "Treat_Type")
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
plot_model(sa.f4,
type = "pred",
terms = c("l.TFT", "Treat_Type"),
axis.title = c("Time from Last Treatment (log transformed)", "Total Count of Pitch Pine"),
title = "Predicted Counts of Pitch Pine Saplings \n >/= 2.5 cm & < 10 cm DBH",
legend.title = "Treatment Type",
line.size = 1,
show.zeroinf = T,
ci.lvl = NA,
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
# plot model two
plot_model(sa.f2,
type = "pred",
terms = c("l.TFT", "Treat_Type"),
axis.title = c("Time from Last Treatment (log transformed)", "Total Count of Pitch Pine"),
title = "Predicted Counts of Pitch Pine Saplings \n >/= 2.5 cm & < 10 cm DBH",
legend.title = "Treatment Type",
line.size = 1,
show.zeroinf = T,
ci.lvl = NA,
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
plot_model(sa.f2,
type = "pred",
terms = c("l.other", "Treat_Type"),
axis.title = c("Other Saplings (log transformed)", "Total Count of Pitch Pine"),
title = "Predicted Counts of Pitch Pine Saplings \n >/= 2.5 cm & < 10 cm DBH",
legend.title = "Treatment Type",
line.size = 1,
show.zeroinf = T,
ci.lvl = NA,
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
sa_lr1 <- glmmTMB(Total_Tally~ l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Mineral + (1|Site/Plot_No) + offset(l.TFT),
data = sa_lrdata,
family = poisson)
AIC(sa_lr1) #209.5
sa1_sr <- simulateResiduals(sa_lr1, n = 1000, plot = TRUE) #does not pass
testZeroInflation(sa1_sr) #yes
sa_lr2 <- glmmTMB(Total_Tally~ l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Mineral + (1|Site/Plot_No) + offset(l.TFT),
data = sa_lrdata,
ziformula = ~1,
family = poisson)
AIC(sa_lr2) #211.5
sa2_sr <- simulateResiduals(sa_lr2, n = 1000, plot = TRUE) #fails
testZeroInflation(sa2_sr) #yes
sa_lr3 <- glmmTMB(Total_Tally~ l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Mineral + (1|Site/Plot_No) + offset(l.TFT),
data = sa_lrdata,
ziformula = ~Region,
family = nbinom2)
AIC(sa_lr3) #215.3
sa3_sr <- simulateResiduals(sa_lr3, n = 1000, plot = TRUE) # fails
testZeroInflation(sa3_sr) #yes
#for zero inflation, tried ~1, Region, Treat Type, l.TFT, 1+Region,
#nbinom won't work
#maybe I should eliminate some elements from the model and take another look via DHARMa and see if that makes a difference